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5 ■ Abstract 



We study the performance of QCD simulations with dynamical Wilson fermions by combining 
the Hybrid Monte Carlo algorithm with parallel tempering on 10^ and 12^ lattices. In order to 
compare tempered with standard simulations, covariance matrices between sub-ensembles have 
D . to be formulated and evaluated using the general properties of autocorrelations of the parallel 



o 



tempering algorithm. We find that rendering the hopping parameter k dynamical does not lead 



^ ' to an essential improvement. We point out possible reasons for this observation and discuss more 

suitable ways of applying parallel tempering to QCD. 

PACS numbers: 11.15.Ha, 12.38. Gc 
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I. INTRODUCTION 



Improving Hybrid Monte Carlo (HMC) simulations of QCD with dynamical fermions 
is a long standing problem. While better decorrelation is, of course, highly desirable for 
all observables, it is of crucial importance for the ones which are sensitive to topological 
sectors. In fact, it has been observed that the r( correlator is definitely dependent on the 
topological charge Q |1|]. Thus it is quite important to look for simulation methods that 
produce realistic Q-distributions. From this point of view the topological charge appears to 
be a good touchstone when looking for improvements by new methods. 

For staggered fermions an insufficient tunneling rate of the topological charge Q has been 
observed 0, |^ . For Wilson fermions the tunneling rate has been claimed to be adequate in 
many cases 0, H. However, since the comparison is somewhat subtle, the reason for this 
could also be that one is not as far in the critical region as with staggered fermions. One 
could fear that simulating closer to the chiral limit, insufficient tunneling could become for 
Wilson fermions as severe as for staggered ones. Indeed, for Wilson fermions on large lattices 
and for large values of k near the chiral limit the distribution of Q is not symmetric even 
after more than 3000 trajectories (see e.g. Figure 1 of Ref. ||^). 

In the method of simulated tempering first proposed in Ref. the inverse temperature 
is made a dynamical variable in the simulations. More generally, any parameter in the 
action can be made dynamical. Let us suppose that, depending on this particular cou- 
pling parameter, the chosen algorithm has a largely different tunneling rate between certain 
metastable states (in configuration space). Augmenting the algorithm with the tempering 
method means that now the system is updated in an enlarged configuration space includ- 
ing the coupling. Instead of overcoming a high barrier at an unfortunate parameter value, 
a detour in parameter space is now opened to be an easier route. This results in better 
decorrelation. 

Considerable improvements have been obtained with dynamical number of the degrees of 
freedom in the Potts- Model [0, with dynamical inverse temperature for spin glass ^ and 
with dynamical monopole coupling in U(l) lattice theory 0. With dynamical mass of stag- 
gered fermions in full QCD ||10| it has been indicated that by tempering a better sampling of 
the configuration space (with respect to topological charge) can be achieved. However, sim- 
ulated tempering requires the determination of a weight function in the generalized action. 
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and an efficient method of estimating it p, ^ turns out crucial for successfully accelerating 
the simulation. 

A major progress was the proposal of the parallel tempering method (PT) |Tl|, [1^ 



m 

which no weight function needs to be determined. This method has allowed large improve- 



ments in the case of spin glass |TT| . In QCD improvements have been reported with staggered 
fermions |]TB||, applying PT to subensembles characterized by different values of the quark 
mass. This has led to the expectation that, analogously in the case of Wilson fermions, 
introducing various values of the hopping parameter k might be the right choice of param- 
eter for applying the idea of PT. In a first study of this problem, simulations of QCD with 



0(a)-improved Wilson fermions [14|, no computational advantage has been found. Because 
only two coupled sub-ensembles, both at relatively small k, had been used, this could not 
be the final answer concerning the potential capabilities of the PT method. In a previous 
work |[T^, with more ensembles and (standard) Wilson fermions on an 8^ lattice, we have 
observed a considerable increase of the transitions between topological sectors. We have 



extended this study to larger lattices (10 and 12 ) in ||T6 



In the present paper, our task will be to compare standard HMC with HMC combined 
with PT in a more elaborate, quantitative manner. In order to really compare algorithms one 
has to relate the computational effort (computer time) to errors of final results (e.g. particle 
masses). In the case of PT the calculation of errors becomes more complicated, because one 
has to take cross correlations between ensembles into account. Cross correlations lead to the 
technical problem of calculating full covariance matrices from (auto-) correlation functions. 
In this paper, we carry out such an analysis for the average plaquette and the topological 
charge. 

At the time when we were doing our autocorrelation analysis SESAM published autocor- 



relation results for all their runs |T^. Qualitatively the SESAM results are very similar to 
ours, i.e. they also do not observe a clear increase of autocorrelation times with increasing 
K (lower quark mass) in the interval [0.1560,0.1575]. One might have expected some mass 
dependence from the positive experience with tempering in the case of staggered fermions. 
In the light of this fact, the final conclusion of our analysis, that the PT method with respect 
to the hopping parameter does not lead to an improved performance of the HMC algorithm, 
is not really astonishing. Nevertheless, formulating the tools for the correlation analysis is 
a major part of this paper which should prove useful for possible future applications of PT. 
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The paper is organized as follows. In Section II we describe the method of PT to be 
applied to the HMC algorithm. Our simulation results are discussed in Section ^TT[ In 
Section we give general properties of covariances and autocorrelations which - in view 
of the moderate statistics available - must be employed to fully exploit the autocorrelation 
data. In Section [V| we present and discuss our results on integrated autocorrelation times. 
The respective results for off-diagonal elements of the covariances are presented and discussed 
in Section In Section [VI1| we compare the efficiency of the two simulation algorithms 



(HMC without and with PT) and study the effect of cross correlations. In Section VIII 



we comment on swap acceptance rates and give a new formula for that rate extending the 



one used in Ref. [14|. We conclude in Section [X and point out potentially more promising 



applications of PT to QCD. 



II. PARALLEL TEMPERING 



In standard Monte Carlo simulations one deals with one parameter set a and generates a 
sequence of field configurations F{s), where s denotes the Monte Carlo time. In our case the 
set a includes the physical parameters (3, k and algorithmic parameters like the leapfrog time 
step and the number of time steps per trajectory. One field configuration F{s) comprises 
the gauge field and the pseudo fermion field. 

In the parallel tempering (PT) approach |jTl], |l^ one updates K field configurations 



with u = 1, . . . , K in the same run. The characteristic feature is that the assignment of the 
parameter sets aj with j = 1, . . . , K to the field configurations F^, changes in the course 
of a tempered simulation. The total configuration at time s thus consists of B{s), Fi{s), 
F2{s),..., Fk{s) where the permutation 




B{s)= (1) 




describes the assignment of the field configurations Fj,^.(s)(s) to the parameter sets aj. For 
short this approach is called PT with K ensembles. 

The update of the F^, occurs by the usual HMC procedure using the parameter sets aj 
as assigned at a given time. The update of B is achieved by swapping pairs according to 
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the Metropolis acceptance condition with probabihty 

-Pswap(^,j) = min (l, e"^-^) , 

/\H = H{a„ + H{aj, - - H{a,, (2) 

where H denotes the Hamihonian of the HMC dynamics for the parameter set aj and the 
field configurations F^,.. The total update of the Monte Carlo algorithm, after which its time 
s increases by one, then consists of the updates of all followed by the full update of B 
which consists of a sequence of attempts to swap pairs. 

Detailed balance for the swapping follows from (||). Ergodicity is obtained by updating 
all Fi, and by swapping pairs in such a way that all permutations (P can be reached. There 
remains still the freedom of choosing the succession of the individual steps. All such choices 
lead to legitimate algorithms, which might differ in efficiency. 

Our choice of steps is such that the updates of all F^ and that of B alternate. Our criterion 
for choosing the succession of swapping pairs in the update of B has been to minimize the 
average time it takes for the association of a field configuration to the parameters to travel 
from the lowest to the largest K-value. This has led us to swap pairs belonging to neighboring 
K-values and to proceed with this from smaller to larger K-values. 

The observables of interest are associated to a specific parameter set aj. We denote them 

as 

= j = l,...,K. (3) 

III. SIMULATION RESULTS 

We have simulated lattice QCD with standard Wilson fermions and one-plaquette action 
for the gauge fields. In the HMC program the standard conjugate gradient inverter with 
even/odd preconditioning was used. The trajectory length was always 1. The time steps 
were adjusted to get acceptance rates of about 70% in the HMC Metropolis step. In all cases 
(standard HMC, tempered runs for all lattice sizes and ensemble sizes) 1000 trajectories were 
generated, with additional 50-100 trajectories for thermalization. 

We have performed tempered runs using 6 and 7 ensembles, all at /3 = 5.6 on 10^ and 12^ 
lattices, as well as standard HMC runs for comparison. Our simulations cover the K-range 
investigated by SESAM {k = 0.156,01565,0.157,0.1575) g. In a large scale PT simulation 

5 



analogous to that by SESAM one would be interested to use all subensembles (separated in 
K,) anyway, such that additional computational effort seems to be affordable. For instance, 
in the run with 6 ensembles we extended the K-range by adding lower values of k, while in 
the run with 7 ensembles we have used a denser spacing of the K-values. Our /t-values are 
listed in Table |. 

The observables determined were the average plaquette P and the topological charge Q. 
The topological charge was measured using its naively discretized plaquette form after doing 
50 cooling steps of Cabibbo-Marinari type. This method gives multiples of a unit charge. 
The value of the unit charge is close to 1 and has to be determined from the measurements. 

Figure |T] shows typical time series of Q obtained for standard HMC and for tempered 
HMC with 6 and with 7 ensembles. One sees that tempering makes Q fluctuate much 
stronger. Such behavior is indicative for the decreasing of correlations between subsequent 
trajectories. The time series on the 12^ lattice exhibits a richer pattern of transitions as 
that on the 10^ lattice, and the width of the topological-charge distribution increases. But 
the rate of fluctuations decreases with increasing k. 



In our previous investigations [15, |T^ we have considered the mean absolute change of 
Q (called mobility in [Q) to account quantitatively for these rates of fluctuations. However, 
this quantity does not provide a quantitative measure of the computational gain obtainable 
with the tempering method in comparison with standard HMC. 

In the present investigation, we therefore base the comparison on the full account of 
the non-diagonal covariance matrix for different observables to be introduced in Eq. (^. 
The covariance matrix will be calculated from general correlation functions which, for an 
observable O and a number N of updates, are defined as 

N N N 

s=l s'=l s"=l 

For j = k they are the usual autocorrelation functions, while for j ^ k they describe cross 
correlations between different ensembles. 

Typical examples of normalized autocorrelation functions Pj{t) = Rjj{t) / Rjj{0) are pre- 
sented in Figure 0. It can be seen that for the tempered runs the decay is considerably 
faster than for the standard ones. Among the tempering runs it is fastest for the run with 
7 ensembles. In the latter case the remarkable fast decay occurring in the interval t G [0, 1] 
should be noticed. 



6 



With the statistics available, the correlation functions decay below the Monte-Carlo noise 
for relatively small t. Although the much faster decay of the functions for tempering is 
apparent, giving numbers for the autocorrelation times and cross correlations is clearly a 
formidable task. In principle, in order to estimate the autocorrelation times within 10 % 
one would need trajectory numbers higher than ours by roughly an order of magnitude. To 
reach a conclusion about the expected improvement we do not attempt such a precision mea- 
surement of the autocorrelation times. Moreover, we want to apply some a priori knowledge 
on the spectral properties in order to exploit our simulation data in the most efficient way 
possible. Therefore, in the following we first have to elaborate on some theoretical issues 
concerning covariances. 



IV. COVARIANCE MATRIX AND MARKOV SPECTRUM 



We obtain the covariance matrix by using the general correlation function and gen- 
eralizing the derivation given in Ref. |TH| for the case j = k 



1 / t 

\ t=i 



The diagonal elements of (|^) are the variances of Oj which are traditionally written in the 
form 

var(0,) = 2r, , (6) 

introducing the integrated autocorrelation times 

^ t=l 

where pj{t) = Rjj{t)/ Rjj{0). 

When evaluating practical simulations the summation in up to — 1 makes no sense 
since Pj{t) is buried in the Monte Carlo noise already for relatively small t. Therefore, it 



has been proposed ||T8|, to sum up only to some smaller value M of t. However, in 
practice that procedure is not stable against the choice of M and neglecting the rest is a 
bad approximation. The proposal to estimate the neglect by an extrapolation based on the 
t values M and M — 1 pO| is still inaccurate in general. A more satisfying procedure is to 
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describe the rest by a fit function based on tlie (reliable) terms of @ for t < M and on the 
general knowledge about the Markov spectrum. This procedure has led to perfect results in 
other applications 

In order to apply the latter strategy also for determining off-diagonal entries in we 
have to look how spectral properties enter the parallel-tempering case. For such con- 
siderations it is convenient to introduce a Hilbert space [|r^, ^ with an inner product 
(0) x) = Xlc /^(^)'^*(^)x(^)' where /i(C) is the equilibrium distribution of the system and C 
denotes the configurations C = {B, Fi, F2, . . . , F^}. Using this notation we can write the 
expectation values (Oj) = Xlc /^(^)^i(^) ~ (I'^i) two-time correlation functions 

{OMOk{t))= Yl /^(C™)0,(CM)iy*(CM;CM)a(CW)= (8) 

as inner products, where 

1^*(CM;CW) = Yl W^(C[°l;CW)l^(CW;C[2])...W^(C[*-i];CW) (9) 
c[ii,...,c['-ii 

is the t-step transition matrix constructed from the one-step transition matrix W{C] C) of 
the Markov process considered. In the spectral representation 

W = Y^ ^rPr (10) 

r>l 

with eigenvalues Xr and projection operators Pr, one has Ai = 1 for the equilibrium eigenvec- 
tor /i(C) and |Ar| < 1 for the other modes. Obviously only Pi contributes to the stationarity 
relation n{C')W {C ; C) = fJ.{C), and Pi{C',C) = fi{C) follows. With this notation one 
can rewrite 

{0,){0k)= [o„P,Ok) . (11) 
Using (H), (0) and (|lT]) we obtain for the general correlation function 

R,,{t) = {O,i0)O,{t)) - {0,){0,) = J2K{o„PrO,) , (12) 

where, due to the subtraction, the term with Ai = 1 cancels out. We thus get the general 
representation 

Rjk(t) = ajfcr A* with \Xr\ < 1 (13) 
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where only the coefficients ajkr depend on the particular pair of observables while the eigen- 
values Xr are universal and characterizing the chosen simulation algorithm. It is important 
to realize that this also holds for the observables of form (H) used in PT. 



V. AUTOCORRELATION RESULTS 



For the numerical evaluation of @ we apply the method explained in Section using the 
fact that after some time only the slowest mode in ( |T3|) survives. Our method is to sum up 
the simulation data only up to some t before the noisy region and determining the rest of the 
sum from a fit assuming that the fit describes the slowest mode well. For tj the rest typically 
amounts maximally to about 25 %. The proper choice of the fit intervals in t (excluding the 
region of fast contributions and the noisy region) was controlled by inspection of the graphs 
and watching the resulting values. Examples of such fits are shown in Figure ^j. 

In view of the moderate statistics available we additionally have made use of the univer- 
sality of the Markov spectrum implying that in (|13|) for given algorithm and lattice size only 
the coefficients ajkr can vary with the observables. We have verified that a collective fit for 
the whole diagonal with a universal slowest mode gives results comparable to individual fits 
(see Table |). That motivated us to perform collective fits with one single mode to diagonal 
and non-diagonal terms. In fact that method greatly helped to get stable fits which will be 



further discussed in section VI 



To obtain errors for the covariances one can generalize the derivations given in Ref. 
T8| for the diagonal case to calculate covariances of covariances from the Rjk{t) data only. 



However, such calculation is impractical with the statistics available. Therefore, we have to 
rely on the comparison of measurements of covariances at different parameter values and on 
consistency checks to get some idea about the size of the errors. 

Comparing the results for Rjj{0) from the different simulation algorithms, which should 
give the same numbers, one sees that this ingredient of a calculation of Tj has errors of 
about 20 %. The exponential autocorrelation time r^^p = — 1/lnA, where A denotes the 
slowest mode, corresponds within a good approximation to the integrated autocorrelation 
time which one would get taking only the slowest mode into account. The presence of faster 
modes then renders tj smaller than Xgxp • We always find consistency with this requirement. 

The integrated autocorrelation times tj obtained in this way are given in Table |. One can 



9 



see that there is httle difference between the results from individual and collective fits. The 
fluctuations of values found for neighboring /t-values indicate relatively large errors. Judging 
from the observed noise levels the errors are expected to be largest for the standard case 
and smallest for tempering with 7 ensembles. Despite these errors two unexpected features 
are clearly visible: i) there is no sizable increase in tj with k and ii) there is nevertheless 
gain in terms of tj when using tempering. 

The lack of a sizable increase in tj with k, which one would have expected for the 
standard runs, firstly indicates that valuing time histories by eye (as of Q in Figure |l|) can 
be misleading. It secondly shows the important fact that the usual precondition of successful 
tempering, connecting regions with considerably different r, is not fulfilled here. 

In the light of this it comes with some surprise that nevertheless gain in terms of Tj is 
observed. The reduction of tj for Q turns out to be larger than for P. 

VI. OFF-DIAGONAL COVARIANCES 

For the discussion of cross correlations in the following Section we need the off-diagonal 
elements of the covariance matrices. As in the diagonal case, the use of the simulation data 
in the sum makes only sense outside the noisy region. One finds that the off-diagonal 
elements of the general correlation functions are decreasing with the distance from the 
diagonal. Therefore, the evaluation of the sum in the off-diagonal case is more difficult 
because beyond some distance \j — k\ the elements are completely indiscernible within the 
noise. See Figure ^ for an illustration. 

In tempering with 7 ensembles we generally find that three off-diagonals can be deter- 
mined while for tempering with 6 ensembles only one. For the off-diagonal elements Rjk{t) 
which clearly show a signal above the noise we generally observe a maximum at t = \j — k\ 
(see Figure ^). We therefore look for a prediction of their functional form. Rewriting (|T2|) as 
Rjk{t) = (^Oj{s), (1 — Pi)W^Ok{s)j their qualitative behavior can be discussed. Obviously 
at each time the observables of type (^ depend only on one field configuration. This form 
of Rjk{t) suggests that for j ^ k a sizable contribution only arises when, under the action 
of the transition matrix by W^Ok{s), a contribution also depending on the field configura- 
tion entering Oj{s) has been generated. For our type of swapping this situation occurs for 
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t ^ \j ~ k\ so that 

{V dikrK- for \j — k\ < t 
for < t < |j - k\ 

should be the approximate behavior. We indeed generally see this behavior within errors in 
our data. 

For the numerical evaluation of (|^) we again apply the method explained in Section 
restricted to a t interval where the respective Rjk{t) signal is sufficiently above the 
noise. To get stable off-diagonal results we use the method of the collective fit, using the 
existence of a universal slowest mode described in Section |V[ Table y shows an example 
of a numerical result for Cjk (remember that Cjk is symmetric). Generally the off-diagonal 
elements obtained, especially the smaller ones, are likely to be overestimated because of 
possible contributions of the noise. 



In the applications to be described in Section |VII| the full covariance matrix Cjk is needed. 
This excludes, for this purpose, the consideration of the PT results in our case of 6 ensembles, 
because we were able to determine elements only in the ffist subdiagonal \j — k\ = 1. They 
have much larger errors than in our PT studies with 7 ensembles. In this case, in contrast, 
we got 3 sub- diagonals. We find that putting the remaining ones equal to zero or using 
various extrapolations in \j — k\ for them makes only little difference in the results. This 
reflects the fact that, although the elements close to the diagonal are not small, there is 
nevertheless a faster decrease farther away from the diagonal. 



VII. COMPARISON OF SIMULATION ALGORITHMS 

At the end of Section [V| we have pointed out that there is no increase of autocorrelation 
times with n. Because of this observation the usual mechanism of tempering - which is to 
provide an easier detour through parameter space for the suppressed transitions - is not 
available here. If this mechanism is working, tempering with several parameter points can 
be advantageous even if one is interested only in the result at one point. This holds, in 
particular, for systems where in the region of interest otherwise almost no transition occurs 
0, ||. Unfortunately our region is not of this type. 

As also observed in Section there is nevertheless a reduction of autocorrelation times. 
The effect of this is, however, not large enough to get generally gain if one is interested in 
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only one point. In fact, dividing the reduction factors of the tj in Table | by the number of 
ensembles it can be seen that at best in the case of Q some gain remains. 

We now turn to the question whether gain remains if one is interested in the results at 
all parameter values. In this case it is necessary to account for cross correlations between 
the ensembles. To be able to do this one has to rely on fits to the data. The respective fit 
method is well know from the treatment of indirect measurements (see e.g. Ref. P3[). For 
proper comparison this method, which leads to improved errors, has to be applied to the 
tempering case as well as to the standard case. 

Final results then are obtained from fits to mean values from individual ensembles. In 
the case of PT the full covariance matrix enters the fit. Although there are difficulties to 
account for the full matrix numerically, we have tried to develop some feeling for its influence 
by making fits to the observables we have measured. 

It is known that (Q) = 0. Therefore our fit ansatz for (Q) is a constant, i.e., the 
fit procedure is a weighted mean using the full covariance matrix. For the plaquette we 
observed that our data are consistent with a linear dependence on k. Therefore we have 
used a linear fit ansatz in k. 

In the following we outline the fit method for the case of the plaquette. The linear fit 
ansatz just mentioned is (P) = xi + X2H, where xi and X2 are the fit parameters. In matrix 
notation we have r] = —Ax with 



A = - 



1 ka- 



X 




(15) 



We also introduce the vector of measured values p = {Pi, . . . , Pk)^ and call the correspond- 
ing covariance matrix Cp. The result of the fit is the minimum of [r] — p)^Cp^{r] — p) which 
lies at x = —{A'^Cp^A)~^A^C~^p. The errors of the result are square roots of the diagonal 
entries of = {A'^C'^A)"^. For fitting a constant to measured values of (Q) corresponding 
formulae apply with X2 = 0. 

Actually we are interested in the variances of the fit-function values. To get them we 
insert x into the fit function, 



r] = —Ax , 



(16) 
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and using (|T6D in the transformation law of covariances we obtain 

Cf^ = ACiA^ . (17) 



The diagonal elements of ( |T7| ) are the variances of interest and the square roots of them the 
improved errors. 

For standard HMC, in these calculations we have used the measurements for the 5 selected 
K values (see Table |) . In the case of tempering the number of measurements is equal to the 
number K of ensembles. 



Table |I| gives data with usual statistical errors and fit results (|T6D with improved errors 
(p!7|). We denote the errors by cq for the usual statistical errors (in column 2), by e/ for 
the improved errors taking into account the full covariance matrix (column 3), and by 
for those obtained only with the diagonal elements of the covariance matrix (column 4), 
respectively. 

The factor {ef/edY describes the influence of cross correlations. We typically find values of 
about 2 to 3 for it. As compared to the reduction factors apparent from Table | this appears 
not large. However, one has to be aware that proper comparison here needs consideration 
of the improved errors in the standard and in the tempering case. 

The relevant computational gain factor for the comparison standard vs. tempering case 
is given by the improved errors and the numbers of ensembles as 

/ standard I tempering \2 nj- / nj- /-|q\ 

\^d ' / ' -'■^ standard/ -'■^ tempering • \ 



In the example in Table |T| this factor appears close to one. However, because of the 
inaccuracy of the standard data the respective fit results are not reliable (giving factors 
from about 0.5 to 4 in other cases). Thus we are not able to give definite numerical results 
for dll). 

A further quantity to be mentioned is the reduction factor Co/ej for the errors. This 
reduction in general is larger for larger cq. Thus the comparison based on the fit procedure 
favors the standard case, which reduces a possible gain. 

Altogether it looks that even with more accurate data it might be difficult to present 
evidence for computational gain in our case, where the autocorrelation times for standard 
HMC do not vary within the K-range considered. 
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VIII. SWAP ACCEPTANCE RATES 



For the effect of enlarging the lattice size on the efficiency of the algorithm the change 
of the swap acceptance rate (Pswap) resulting from appears to be relevant. In Ref. fT^ 



agreement has been reported between the behavior observed in PT with the expression 
erfc ^Iy^ (AH)^ derived in Ref. [24|. The derivation given there, however, relies on the 
area-preserving property of the HMC algorithm, implying (exp(— AiJ)) = 1, which does not 
hold in the case of swapping. 

To get the relation appropriate for swapping, we again neglect higher cumulants in the 
cumulant expansion {exp{—AH)) = exp ^ — {AH) + ^{{AH — {AH))'^)) =1= . . . j . However, 
in contrast to Ref. [^, we put (exp(— Aif)) = exp{—6) with an unknown 6. By convexity 
of the exponential function one finds that (AH) > 6 holds. The relation between mean and 



width of the gaussian used in Ref. p4| then generalizes to 

{AH) = ^a' + 6 , a' = {{AH~{AH)f). (19) 
The evaluation of the integral ^^7=^ dx min(l, exp(— x)) exp ^ — ^'^^^^^^^ j gives then 

where u = a/ {AH) — 6. 

Within errors our values of {AH) turn out to scale with the volume (being roughly 
1.4 and 2.8 for 6 ensembles and 0.35 and 0.7 for 7 ensembles for 10^ and 12^, respectively). 
While the values of {exp{—AH)) for the 10'^ lattice within errors conform with 1, on the 12^ 
lattice they deviate substantially from this (increasing from 0.6 to 1.4). 

Our data agree with (^) using 6 as found from our simulations. In the cases where 
we find that 5 = is not true, despite of this using 5 = in the acceptance formulae 
within errors gives still consistency. However, for larger lattices, where further increase of 
the deviations of {exp{—AH)) from 1 is to be expected, this might be no longer so. The 
indicated deviations tend to improve the situation on larger lattices. Since the behavior of 
6 is not known quantitatively, detailed predictions on {AH) at present are not possible. 
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IX. CONCLUSIONS 



In this paper we have compared PT with standard HMC quantitatively on 10^ and 12^ 
lattices. We have described the steps of such an analysis and carried them out for the 
average plaquette and the topological charge. In a quantitative analysis one has to look at 
the size of errors and therefore in principle also at errors of the errors. While the first is 
at the limits of feasibility, the latter is definitely beyond our statistics. For that purpose 
we have demonstrated how the cross-correlation functions between the sub-ensembles have 
to be taken into account, using their general, algorithm independent properties. This has 
allowed us to make consistency checks, and we found consistent behavior of the results. We 
believe that this part of our analysis is of general interest. 

The choice of the /t-range in this paper was guided by recent large-scale QCD simulations 
with dynamical Wilson fermions like SESAM [^] and triggered by the suppression of tun- 
neling rates for the topological charge observed at the largest K-values available. However, 
this range might be still too far from the chiral limit with the consequence of no dramatic 
change of the autocorrelation times for the standard Hybrid Monte Carlo method. This 
would explain why PT in our case, i.e. for the k range covered by the work of SESAM, did 
not provide a considerable computational gain. 

At stronger coupling we know that approaching the chiral limit we arrive at a second order 
transition into a phase with broken combined parity-flavor symmetry (so-called Aoki-phase) 
| 25| . At this transition we expect a critical slowing down and therefore strongly increasing 



autocorrelations. This is the reason, why we nevertheless believe that PT should become 
efficient if one is sufficiently close to the chiral limit. An even more promising scenario might 
be the application of the PT method towards the continuum limit, in particular along lines of 
constant physics in the (3 — K-plane. Already in pure SU (3)-Yang-Mills theory the tunneling 
rate between different topological sectors becomes strongly suppressed with increasing jS. 
The tools developed in the present paper will be very useful for such future applications. 
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TABLE I: Integrated autocorrelation times Tint,p for plaquette and Tint,Q for topofogical charge 
from individual fits (see Section^. For comparison results from collective fits are given in brackets. 
For each case 1000 trajectories were generated with trajectory length 1. 

standard HMC tempered HMC 

6 ensembles 7 ensembles 











= 0.0005 


Ak = 


0.00025 




10^ 


12^ 


10^ 


12^ 


10^ 


12^ 


n 








2rint,p 






0.15500 


14.0 


16.7 


5.9 [6.7] 


8.8 [9.5] 






0.15550 






6.3 [6.9] 


10.2 [8.6] 






0.15600 


10.4 


13.3 


3.7 [4.0] 


5.5 [6.2] 


4.2 [4.6] 


4.6 [5.1] 


0.15625 










2.8 [2.6] 


2.8 [3.4] 


0.15650 


9.3 


12.9 


5.2 [5.6] 


6.1 [6.9] 


2.9 [2.5] 


3.4 [3.4] 


0.15675 










2.7 [2.4] 


4.5 [3.7] 


0.15700 


8.6 


14.2 


5.8 [5.6] 


7.8 [8.1] 


2.9 [2.4] 


6.1 [4.8] 


0.15725 










4.1 [2.9] 


5.2 [4.9] 


0.15750 


9.0 


8.1 


8.4 [8.5] 


10.9 [10.8] 


3.4 [3.9] 


8.2 [7.9] 










2rint,Q 






0.15500 


42 


22 


15 [17] 


17 [23] 






0.15550 






13 [11] 


11 [12] 






0.15600 


37 


74 


7 [7] 


20 [19] 


16.6 [16.5] 


4.8 [7.7] 


0.15625 










9.7 [9.7] 


5.2 [6.3] 


0.15650 


41 


48 


8 [7] 


16 [17] 


5.7 [6.1] 


5.8 [6.5] 


0.15675 










3.3 [3.1] 


6.4 [6.8] 


0.15700 


45 


38 


16 [12] 


8 [10] 


2.3 [2.5] 


6.3 [6.5] 


0.15725 










5.3 [4.9] 


5.3 [5.5] 


0.15750 


46 


14 


6 [8] 


35 [26] 


6.8 [7.6] 


5.5 [7.0] 
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TABLE II: Covariancc matrix of Q for tempering with 7 ensembles on 12^ lattice. 



J 


Kj 


W ■ Rjj 


103 . i?,- 






1 


0.15600 


12.68 


10.34 


6.87 


3.14 


2 


0.15625 


9.89 


8.28 


5.77 


2.04 


3 


0.15650 


9.90 


8.44 


4.42 


1.04 


4 


0.15675 


9.27 


6.30 


3.22 


1.23 


5 


0.15700 


7.02 


4.46 


2.39 




6 


0.15725 


4.34 


2.82 






7 


0.15750 


4.26 









TABLE III: Data and fit results for P from standard runs and from tempering with 7 ensembles 
on the 12^ lattice. 



K 


data 


fit with full 
covariance matrix 


fit with diagonal 
covariance matrix 


Standard KMC 


0.15500 


0.431825 (135) 




0.432002 (112) 


0.15600 


0.429872 (136) 




0.429968 (66) 


0.15650 


0.429343 (132) 




0.428952 (55) 


0.15700 


0.428309 (138) 




0.427935 (59) 


0.15750 


0.426695 (91) 




0.426919 (77) 


tempering (7 ensembles) 


0.15600 


0.430024 (79) 


0.429973 (69) 


0.430064 (50) 


0.15625 


0.429494 (64) 


0.429516 (56) 


0.429578 (38) 


0.15650 


0.429138 (63) 


0.429060 (48) 


0.429093 (30) 


0.15675 


0.428749 (67) 


0.428603 (48) 


0.428607 (29) 


0.15700 


0.428166 (81) 


0.428146 (56) 


0.428121 (35) 


0.15725 


0.427560 (82) 


0.427690 (68) 


0.427636 (46) 


0.15750 


0.427034 (102) 


0.427233 (84) 


0.427150 (60) 
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standard HMC 



Tempered HMC 

6 ensembles 
Ak = 0.0005 



Tempered HMC 

7 ensembles 
Ak = 0.00025 



K = 0.1565 _ 



K = 0.1565 



,11 1 .11' 



II 



K = 0.1565 






1 ' 


= 0.1575 









K = 0.157 



iiY-'--vp|ftl|j 



Q - 



K = 0.1575 



b4 



600 traj. 



1000 



traj. 




K = 0.1575 



1000 



500 traj. 1000 



FIG. 1: Time series of Q for standard and tempered HMC on 12^ lattice at /? = 5.6 (for part of 
K- values only, see Table | for full list). 
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FIG. 2: Normalized autocorrelation functions for Q for standard HMC (left) and tempering with 
6 ensembles (center) and 7 ensembles (right) on a 12^ lattice for (3 = 5.6, k = 0.1565. The errors 
indicated are the purely statistical ones. The lines represent fits to the subset of data points with 
full symbols. 




1.5 



^ 



FIG. 3: Correlation functions for Q for tempering with 7 ensembles on a 12^ lattice for /3 = 5.6. 
Shown are the autocorrelation function at k = 0.1565 (left) and the cross-correlation functions 
for K = 0.1565,0.15675 (center) and k = 0.1565,0.157 (right). The errors indicated are purely 

statistical ones. The lines represent one combined fit as explained in the text. Full data symbols 
indicate the fit interval. 
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